VISCOUS VORTEX FLOWS 


N88 - 14934^ 

~C> 2 ^ 

// 7<?3 2 — 


R. P. Weston, J. P. Chamberlain, 

C. H. Liu, and Peter-Michael Hartwich* 
NASA Langley Research Center 
Hampton, Virginia 


*NRC - NASA Resident Research Associate 



tODOCTIOH 


Several computational studies are currently being pursued that focus on 
various aspects of representing the entire lifetime of the viscous trailing 
vortex wakes generated by an aircraft. The computational regions are 
indicated in the figure below, with the vicinity around the vortex-generating 
aircraft designated as region A. The formation and subsequent near-wing 
development of the leading-edge vortices formed by a delta wing are being 
calculated at modest Reynolds numbers using a three-dimensional, time- 
dependent Navier-Stokes code. The calculations exhibit realistic vortex 
characteristics including behavior similar to vortex bursting at high angles 
of attack. Another computer code has been developed to focus on the roll-up, 
the trajectory, and the mutual interaction of the trailing vortices further 
downstream from the wing (region B) using a two-dimensional, time-dependent 
Navier-Stokes algorithm. This code has also been used to study the 
modification of the vortex behavior due to ground-proximity effects and the 
enhanced vortex decay induced by atmospheric temperature gradients. To 
investigate the effect of a cross-wind ground shear flow on the drift and 
decay of the far-field trailing vortices (region C), yet another code has been 
developed that employs Euler equations along with matched asymptotic solutions 
for the decaying vortex filaments. And finally, to simulate the conditions 
far downstream after the onset of the Crow instability in the vortex wake 
(region D), a full three-dimensional, time-dependent Navier-Stokes code has 
been developed to study the behavior of interacting vortex rings. 
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3-D NAVIER-STOKES CALCOLAffQNS 

The formation and the subsequent near-wing development of the leading- 
edge vortices formed by a delta wing are simulated by time-accurate finite- 
difference solutions to the Navier-Stokes equations for three-dimensional, 
incompressible flows with moderate Reynolds numbers. This particular aspect 
of an extensive numerical study of vortex flow fields generated by aircraft 
was chosen because a broad experimental data base allows thorough evaluation 
of the numerical results. The discretized momentum equations are integrated 
by an Euler-explicit, time-marching procedure that is first and second-order 
accurate in time and space, respectively. The pressure field is computed 
simultaneously by a simplified Jacobi method applied to the Poisson equation 
for the pressure. Truly transient solutions are achieved by repeating the 
iterative procedure as with time level until the flow is source free. 


• Body-fitted coordinates 
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COMPARISON WITH EXPERIMENT AT a - 20.5° 


This figure shows a comparison of streamline patterns that were 
determined from the velocity fields as steady-state solutions to the Navier- 
Stokes equations (Re = 250 and Re = 1500) and to the Euler 

equations (Re + «>) , and also those that were photographed in a water tunnel 
(Re = 900). In all four cases the lateral deviations of the vortex 
streamlines are quite similar. The circumferential velocities of the 
experimentally observed vortices are judged to be greater than those in the 
computed flow fields since the number of coils per axial length in the 
photograph is higher than in the computations. The values for < give an 
estimation of the additional numerical viscosity which has to be added to 
numerically stabilize the finite-difference solutions of the Navier-Stokes 
equations (Re = 1500) and of the Euler equations (Re ^ °°) . The quantity k can 
be interpreted as the reciprocal of a Reynolds number, giving the effective 
Reynolds number during the computations as <^> effective - ^nominal + K * 






Re = 900 experiment 




Re=900 experiment 
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VELOCITY FIELD FOR 78.6% CHORD, a = 20.5° 

This figure shows the projections of the velocities in planes which 
intersect the delta wing perpendicularly at 78.6% of the root chord. The 
aspect ratio is AR = 1 and the angle of attack is 20.5 degrees. The velocity 
fields represent steady-state solutions to the Navier-Stokes equations (Re = 
250 and Re = 1500) and to the Euler equations (Re -* <») • For Re = 250, no 
additional numerical damping was necessary, whereas numerically stable 
solutions for Re = 1500 and Re + °° were obtained only after additional 
damping, indicated by k, has been introduced in the discretized momentum 
equations. The vortex for Re = 250 lies higher and closer to the plane of 
symmetry than in the other two cases. The velocity distributions for Re = 

1500 and Re * «> look quite similar if a small region in the neighborhood of 
the wing is neglected where the different boundary conditions for viscous and 
inviscid flows affect the solution. The small size of this region suggests 
that for the computation of vortex flows around sharp-edged delta wings with 
Re > 0(10°) the friction forces may be neglected. To provide further evidence 
for this conclusion, additional investigations are in progress. 



Re=250 Re = 1500 k=1.8-10" 4 Re^oo k =74-10' 4 
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TWO-DIMENSIONAL INITIAL VALUE PROBLEM 


The governing equations are expressed in terms of the variables of stream 
function ( \p) and vorticity ( a>) . The equation of continuity can then be cast 
as a Poisson equation and the curl of the momentum equation becomes the 
vorticity transport equation, where t is time, v is the kinematic viscosity, 
and v and w are the velocity components in the y and z directions, 
respectively. The subscripts denote partial derivatives. The initial 
vorticity distribution provides the initial condition and the boundary 
condition for the unbounded flow is based on the fact that the vorticity 
distribution is confined to one area and decays exponentially with distance 
from that area. The exact solution for the Poisson equation is given by the 
Biot-Savart law. However, the computation of the boundary conditions using 
this law directly is extremely expensive. A considerable time savings can be 
achieved by expanding the Biot-Savart integral in a power series in terms of 
vorticity. The coefficients of the terms in the power series are moments of 
the vorticity distribution. Previous studies have shown that the first moments 
and several linear combinations of higher moments are time invariant. These 
results can be used to determine the boundary condition for the numerical 
calculations and are also used in checking the accuracy of the numerical 
solutions as they proceed. 

I ncompress i b I e , Nav i er-S tokes Eqns . 

Continuity: V*V = 0 => A ^ = -« (Poisson eqn) 

Vor t i c i ty 

Transport : + v w <y 2 = v A w 

I.C. w ( x , 0 ) = w 0 (x) 

B.C. w(x,t) -> 0 exponentially as x oo 

Boundary values of f are evaluated using 
a moment expansion of the Biot-Savart law 
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FLOWCHART 


The basic algorithm performed by the program is outlined below: 

1 . The volume integrals are evaluated using Simpson integration over the 
computational domain. The time varying integrals are evaluated at each 
time step, while the time-invariant integrals are evaluated initially and 
then only periodically to monitor the accuracy of the solution. 

2. Ihe results from step one are used to obtain the expansion coefficients 
which are in turn used to determine values of at the boundary. 

3. A fast Poisson solver is used to determine the values of in the interior 
of the computational domain. The solver currently in use is a direct 
method developed by the National Center for Atmospheric Research (NCAR) . 
The method uses a finite-difference formulation and is second-order 
accurate in the spatial directions. 

4. The velocity field is obtained by using second-order centered differences 

to obtain the curl of The velocity values are written over the stream 

function values to minimize the required computer storage. 

5. The vorticity field is advanced in time by using a finite-difference 
representation of the vorticity transport equation. The program currently 
uses the Dufort-Frankel method (ref. 1) for solving the vorticity 
transport equation; this explicit method is accurate to 

0[(At) 2 , ( Ax. ) 2 , v( At/Ax i ) 2 ; i = 1, 2, 3]. Explicit methods for solving 
the vorticity transport equation appear to be more appropriate than 
implicit methods in this case. 
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EVOLUTION OF VORTICITY FOR SPANLOAD OF TRANSPORT IN LANDING 

CONFIGURATION 

This figure displays the resultant contours of vorticity for successive 
planes downstream of the wing. The calculations were performed at a Reynolds 
number (T/v) of 20,000. The results are shown for only the right half-plane 
since the flow is assumed to exhibit mirror symmetry about the line y = 0, 
Since the computational grid follows the vortices, the projection of the wing 
trailing edge appears as the horizontal line that rises in successive 
planes . 

The example is for a spanload distribution like that obtained with a 
transport aircraft using flaps on landing or takeoff. The resulting initial 
vorticity distribution has been simplified by distributing the vorticity along 
a horizontal line, unlike the more complicated positions expected in 
reality. Note the negative vorticity values associated with the reduction in 
lift around the fuselage-wing junction. 

Successive frames illustrate the evolution of the vortex system in time 
and demonstrate the ultimate merging of the vortices from the wing tip and 
flap tip. A movie generated from these calculations has been useful in 
observing the progress of these vortex interactions. 
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INTERACTION OF DECAYING TRAILING VORTICES 
IN SPANN IS E SHEAR FLOW 

The drift of trailing vortices in a crosswind near to the ground is 
modeled by an unsteady, two-dimensional, rotational flow field with a 
concentration of large vorticity in spots having a finite but small effective 
size and finite total strength. The problem is analyzed by a combination of 
the method of matched asymptotic analyses for the decay of the vortical spots 
and the Euler solution for the unsteady rotational flow. Using the method of 
averaging, a numerical scheme is developed in which the grid size and time 
step depend only on the length and velocity scales of the background flow and 
are independent of the effective core size of a vortical spot. The core size 
can be much smaller than the grid size while the peak velocity in the core is 
inversely proportional to the spot size. Numerical results are presented to 
demonstrate the strong interaction between the trajectories of the vortical 
spots and the redistribution of vorticity in the background flow field (ref. 
2 ). 



COORDINATE TRANSFORMATION 


The problem of a steady far-wake vortex can be simplified by reducing the 

problem to an equivalent unsteady two-dimensional problem in a plane normal to 

the flight direction* This simplification is performed by changing the 

coordinate z to a time variable t using z=z*-W a> t, where z is stationary and z* 

points in the downstream direction and moves with the airplane at 

velocity W (see figure). This assumption ignores the streamwise curvature of 
00 

the trailing vortex filaments, their initiation at the trailing edge, and the 
variation of the velocity parallel to the z axis. Mathematically, we assume 
that d/dz << d/dx and d/dy with x,y as the span and vertical direction, 
respectively. In the x-y plane at a station z, the trailing vortices are 
represented by "vortex spots" of small effective size inside of which there is 
a strong vorticity distribution with finite total strength. 

This model is employed to study the drift and decay of farfield trailing 
vortices (vortical spots) in a cross wind (a spanwise shear flow) near the 
ground • 



I 
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EFFECT OF VORTEX STRENGTH ON TRAJECTORY 

To gain a qualitative understanding of the interaction of vortical spots 

with a background shear flow, we study first the case of a single vortical 

spot. This figure shows the trajectories of a single concentrated decaying 

vortical spot of various strengths submerged in a background shear flow. The 

initial vertical position of the spot is at y=1 and the initial background 

shear flow is chosen as U^(y)=1-e 1 • The results show that the vortical spots 

with positive circulation drift downstream and upward while the vortical spots 

with negative circulation drift downstream and downward but eventually turn 

backward (t>t*). This phenomenon is more pronounced as the strength of the 

vortical spot increases. To explain this phenomenon we consider the case of a 

single vortical spot with T>0. The disturbed flow moves downward behind the 

spot, x<X, and upward ahead of the spot, x>X. For an initial background 

vorticity with a)^'(y) >0, the disturbed flow increases the vorticity behind 

the spot and decreases the vorticity ahead of it, i.e., £>0 for x<X 

and £<0 for x>X. The background vorticity variation £ in turn induces an 

upward motion of the vortical spot for T>0. From similar arguments, we can 

explain that the background vorticity variation £ will induce a downward 

motion of the vortical spot with T<0. The reason that a vortical spot of 

negative strength turns around and drifts upstream as it gets closer and 

closer to the ground can be attributed to the decrease of the contribution of 

the background shear flow to the forward velocity of the spot and to the 

increase of the induced velocity by the image of the vortical spot with 

respect to the ground, y=0. It should be pointed out here once more that the 

vortical spot will drift horizontally when the background shear flow is either 

a uniform flow (a) =0) or a constant shear flow ( w ^constant) and there will be 
0 ~ u 

no change in the background flow, £=0. 

Same initial height 

r 
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VORTICAL PAIRS SIMULATING TRAILING VORTICES 
IN A CROSSWIND 

The trailing vortex wakes far downstream of an aircraft are modeled by a 
simple vortex pair whose vorticity distributions are concentrated and are 
centered at ( +X^ (0) , Y^fO) ) with strength + r and effective vortical size 6(0). 
The goal of the following numerical examples is to simulate the interaction of 
the decaying trailing vortical pairs subjected to a cross-flow (spanwise) 
ground shear* The background shear flow used in the examples is an 
exponential profile. 

In order to find out when we have to use the shear layer solution, we 
studied the trajectories of a pair of vortices in a shear layer for different 
initial vortex heights. The Y R min approaches an asymptotic value of 2.6 when 
Y k (0) is greater than 7.0. This means that when the vortex spots are above 
y=7, they are far above the shear layer and the interaction with the shear 
layer is negligible. The corresponding trajectories (in real spanwise 
positions) of the vortical pair, starting at different heights Y^( 0)=1 , 2, 3, 4, 5 
in the shear layer, are displayed in this figure to show that the trajectories 
of the vortical spots are sensitive to the starting height, i.e., the altitude 
of the airplane relative to the thickness of the shear layer. 


EFFECT OF DIFFERENT INITIAL POSITIONS 
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NAVIER-STOKES CALCULATIONS FOR UNSTEADY, 

THREE-DIMENSIONAL VORTEX-DOMINATED FLOWS 

A finite-difference Navier-Stokes code has been developed for calculating 
unsteady, three-dimensional, vortex-dominated flows in unbounded fluid 
domains. The algorithm uses an improved boundary condition specification 
which allows the unbounded nature of the physical problem to be represented on 
a finite computational domain. This boundary condition specification permits 
the efficient computation of flows due to closed ring-like vortical tubes or 
structures. These structures are important elements in fluid flows such as 
free jets, atmospheric turbulence, and the far-field wakes of aircraft, and 
studies of their interaction may aid in an understanding of complex vortical 
fluid flows. 

The primary variables used in the computations are the vorticity and the 
vector velocity potential, which is defined as a divergence-free vector field 
whose curl yields the velocity. This definition of the vector velocity 
potential automatically satisfies the incompressible continuity equation, and 
for unbounded flow relates the vector potential directly to the vorticity 
through an integral relationship (the vector Poisson integral). In theory, 
this integral relationship could be used to yield the velocity directly from 
the vorticity, but the numerical evaluation of the Poisson integral throughout 
the computational domain is very time consuming. The current algorithm avoids 
the expense of directly computing the integral by approximating the integral 
values at the domain boundary with a series representation, and then using 
these values along with a fast Poisson solver in the domain interior. This 
technique yields a fast, accurate solution for the velocity field of the 
unbounded-domain physical problem with a finite computational domain. 

Incompressible, Laminar Flow 


Continuity: VV = 0 => A A = -&T (Poisson eqn) 

Vor t i c i ty 

Transport : « t = V x (V x w) + v A u 


I .C. u(x ,0) = w 0 (xj 


B.C. «(x,t) 0 exponentially as 


oo 


Boundary values of A 
a moment expansion of 


are evaluated using 
the Biot-Savart law 
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BOUNDARY CONDITION ACCURACY 


The accuracy and efficiency of two approximate vector-potential boundary 
condition methods ware evaluated by comparing them with an exact solution that 
used Poisson integral evaluations to generate the boundary conditions. The 
two approximate methods were the truncated series method mentioned in the 
previous figure and a method that enforces zero normal velocity at the 
boundary by setting A to zero at the boundary. 

The vorticity distribution chosen for the boundary condition checks is 
shown in the upper right figure as a vorticity-magnitude isosurface. This 
distribution represents two Gaussian-core vortex rings with equal strengths, 
radii, and core diameters whose axes of symmetry lie in the x-y plane and 
cross the x-axis at an angle of ±22.5°. The centers of the rings lie on the 
y-axis at ±1.5 units, where a unit of length is the toroidal radius of each 
ring. The effective core radius (the radius at which the vorticity magnitude 
has fallen to 1 /e of the maximum vorticity magnitude) is 0.5. The computa- 
tional domain is a cube with edges of length 8 centered about the origin. 

This domain size is the practical minimum cubical size that will still enclose 
both vortical rings and thus represents a "worst-case" test condition. 

The three boundary condition figures show contour plots of the magnitude 
of it in the x-y plane of the domain for each of the three boundary condition 
methods. It is evident from the figures that the finite solution domain 
affects both the series method and zero normal-velocity method solutions, 
although the series solution is not affected to nearly as great an extent. 

For slightly larger domains the series method solution closely approximates 
the Poisson method solution, whereas the global character of the zero normal- 
velocity solution is altered by the closed boundary. 



POISSON INTEGRAL B-C-'S 



VORTICITY DISTRIBUTION 



SERIES B-C.'S 



f*l = 0 B.C.'S 
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MERGER OF VORTEX RING PAIR 


This figure shows the oblique collision, merging, and deformation of two 
circular, Gaussian— core vortex rings. The first two plots show representative 
isosurfaces of the initial vorticity and vector velocity potential magnitudes, 
respectively, and the last four plots show the time development of the vector 
potential isosurface. The intersections on the surfaces represent grid point 
locations, and the cube represents the computational domain boundary. The 
rings collide and merge to form a single distorted oblong ring, which 
continues to deform until the major and minor axes have switched their 
orientation. This behavior has been observed experimentally, and work is 
currently in progress to further correlate the code results with experiment. 

It is hoped that this algorithm will be useful in understanding and analyzing 
the physics of vortex-dominated flows. 


(VECTOR POTENTIAL STREAM SURFACES) 



t = 20.2 t = 40-9 t = 60.1 
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